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Inelastic tunneling through magnetically anisotropic molecules is studied theoretically in the presence of a 
strong magnetic field. Since the molecular orientation is not well controlled in tunneling experiments, we con- 
sider arbitrary angles between easy axis and field. This destroys all conservation laws except that of charge, 
leading to a rich fine structure in the differential conductance. Besides single molecules we also study monolay- 
ers of molecules with either aligned or random easy axes. We show that detailed information on the molecular 
transitions and orientations can be obtained from the differential conductance for varying magnetic field. For 
random easy axes, averaging over orientations leads to van Hove singularities in the differential conductance. 
Rate equations in the sequential-tunneling approximation are employed. An efficient approximation for their 
solution for complex molecules is presented. The results are applied to Mni2-based magnetic molecules. 

PACS numbers: 73.63.-b. 75.50.Xx, 85.65.+h, 73.23.Hlc 



I. INTRODUCTION 

Electronic transport through single magnetic molecules has 
recently attracted a lot of interest, 

partly motivated by the hope for applications in molecu- 
lar electronicS j'^-'^'^^ in particular for memory devices and 
quantum computation. On the other hand, the systems also 
show fascinating properties of fundamental interest. Among 
effects observed or predicted for transport through mag- 
netic molecules are the Kondo effectr*^ negative differential 
conductance^^ also at room temperature^S large spin accu- 
mulation in the leads controlled by the initial spin state of the 
molecule,— and large current-induced magnetization changes 
in monolayersji^ 

Molecules with magnetic anisotropy are observed?'^ or 
predicte dVi^-'°i^^ to show particularly rich physics. Magnetic 
molecules based on a Mni20i2 core with organic ligands 
(henceforce Mni2) have a large spin S = 10 in the neutral 
state and strong easy-axis anisotropyi ^'^'^-^-^^ The anisotropy 
leads to an energy barrier between states with large positive 
or negative spin components along the easy (z) axis. Mni2 
has been studied intensively with regard to quantum tunneling 
through this barrieri22i2^ Recently, charge transport through 
Mni2 has been studiedi^*^^ 

Romeike et al^ investigate the effect of quantum tunnel- 
ing on transport through Mni2. They consider small non- 
uniaxial second-order anisotropy terms (5^)^ — (S'^)^ = 
(5^)^/2 + {S~Y /2 and higher-order anisotropics, see also 
Ref. S The non-uniaxial terms are not present for the isolated 
Mni2 molecule^i but are due to the lowering of symmetry by 
the leads. These non-uniaxial terms destroy the conservation 
of the z component of the molecular spin, S^^^. The molecu- 
lar energy eigenstates are then not simultaneous eigenstates of 
S^Qj. Since the non-uniaxial perturbation is small, one can still 
work in a basis of S^^^ eigenstates, but incurs weak tunneling 
transitions between them, which are at the basis of Refs.|a|7|. 
Importantly for the present discussion, they also study Mni2 
in a magnetic field aligned with the easy axis. 

The present paper starts from the realization that there is a 
much stronger symmetry-breaking effect if an external mag- 



netic field is applied. In present-day break-junction and elec- 
tromigration experiments, the orientation of the molecule rel- 
ative to the laboratory frame is not well-controlled. How- 
ever, if the angle 9 between the easy axis and the magnetic 
field is not small and the Zeeman and anisotropy energies are 
comparable — a case that can be realized for Mni2 — there is 
no small parameter. The molecular eigenstates are not ap- 
proximate eigenstates of the spin component along any axis. 
This is the situation studied in the present paper. 

In Sec. the model Hamiltonian is introduced and 
the calculation of the differential conductance is discussed. 
We employ a density-matrix formalism, which allows to 
treat the strong electronic interaction on the molecule 
exactly and is not restricted to low bias voltages, but 
treats the tunneling between molecule and leads as a 
weak perturbation i"^-^-^'^°-'^'^^'^"-^^'^^'^^-^"^-^^'^^ An approxima- 
tion scheme for solving the resulting rate equations is in- 
troduced, which works excellently at low temperatures and 
leads to a large acceleration of the numerics. Results are 
presented and discussed in Sec. |III1 starting with Coulomb- 
diamond plots of differential conductance g vs. gate and bias 
voltages. Then g is shown for varying magnetic field. These 
results are also relevant for situations without a gate such as 
the STM geometry or monolayers of molecules with parallel 
easy axes. Finally, the differential conductance for monolay- 
ers of molecules with random easy-axis orientations is calcu- 
lated. The main results are summarized in Sec. lIVI 

II. MODEL AND THEORY 

A magnetic molecule coupled to two conducting leads L 
(left) and R (right) is described by the Hamiltonian H = 

Hnio\ + H]^ + H-R + -ffhyb- Here, ffmoi is the Hamiltonian 
of the isolated molecule, 

a 

- (if2 + K ^ 4c„) {S^f - H • (s + S) - Js • S, (1) 
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where creates an electron in the lowest unoccupied molec- 
ular orbital (LUMO), s = X^o-o-' "^o- W aa' 1'^ <^a' is the spin 
operator of electrons in the LUMO, and S is the operator of 
the local spin. We also define the total spin Stot = s + S. 
eo is the onsite energy of electrons in the LUMO, which can 
be shifted by applying a gate voltage Vg, U is the repulsion 
between two electrons in the LUMO, and J is the exchange 
interaction between the electron spin and the local spin. A de- 
pendence of the anisotropy energy K2 + kit, on the electron 
number n is taken into account.— 

The g-factors for the local and electron spins are assumed to 
be equal and are absorbed into the magnetic field H, together 
with the Bohr magneton jiB- We choose the easy axis of the 
molecule as the z axis in spin space. Due to the rotational 
symmetry of iJmoh the magnetic field can be assumed to lie 
in the xz plane. We neglect the small non-uniaxial anisotropy 
due to the symmetry breaking by the leads and small higher- 
order anisotropies in order to concentrate on the large effect 
due to the interplay of anisotropy and Zeeman terms. 

The Hamiltonian explicitly takes the electron in the LUMO 
into account, which is non-degenerate for Mni2i^ Since 
the next higher energy orbital (LUMO+1) lies about 8 meV 
above the LUMO,^^ this and higher orbitals are expected to af- 
fect the results only weakly at the low bias voltages discussed 
here. Excitations to one of the higher-energy orbitals are pre- 
sumably responsible for the 14 mV peak observed in Ref.jsl 

We assume that the molecule is symmetrically coupled 
to the leads a = L, R by hybridization terms TJhyb = 

J2aWai^^ly\icr^'^ + l^'C-)' where aly^^ creates an electron in 
lead a. The leads are described by non-interacting electrons, 

TIT _ , „t „ 4.9.10.15 

There are important consequences of the presence of the 
transverse field H^- Since this term does not commute with 
the anisotropy term, 5*^^^. is not conserved and the molecular 
eigenstates (of iJmoi) are not simultaneous eigenstates of S^^^.. 
On the other hand, if the magnetic field is aligned with the 
easy axis^iiSii^ the eigenvalue m of S^^^ is a good quantum 
number and electron tunneling is governed by selection rules 
Am = ±1/2 for sequential tunneling. These selection rules 
do not apply to our case. The only selection rule still vaUd 
stems from the conservation of charge and states that the elec- 
tron number changes by ±1 for sequential tunneling. Apart 
from this, all transitions are allowed. For local spin quan- 
tum number S, there are 25* + 1 molecular states with n = 
electrons in the LUMO and 2(25* + 1) states with n = 1 elec- 
trons. Thus there are 2(25* + 1)^ transitions between states 
with n = and n = 1, leading to a much more complex 
differential conductance g than for the previous case i^i'^i'^ 

The same model can also be used to describe a monolayer 
of magnetic molecules sandwiched between conducting elec- 
trodes, if the tunneling rate between electrodes and molecule 
is sufficiently small to justify the perturbative treatment. This 
can in principle be achieved by spacer groups in the molecules 
or by ultrathin oxide layers. If the interactions between the 
molecules can be neglected, they conduct electrons in paral- 
lel, and the differential conductance per molecule is just the 
ensemble average we are calculating in any case. For a single 
molecule, we have to re-interpret the ensemble average as a 



time average over time scales long compared to the character- 
istic tunneling time. 

These remarks only hold for identical molecules, which re- 
quires all molecules to have their easy axes aligned in par- 
allel. This is a reasonable assumption for molecules that 
are deposited or assembled on the substrate with a preferred 
orientation^^ This has been demonstrated for mixed Mni2 
complexes containing different ligands.=^9 For molecules of 
approximately spherical shape one rather expects the orien- 
tation of the easy axis to be random. To study this case, we 
average the differential conductance over all possible orienta- 
tions. Conversely, measurements of the differential conduc- 
tance can be employed to determine the degree of alignment. 

The absence of constants of motion other than particle num- 
ber requires to diagonalize i7i„oi numerically. The eigenstates 
typically contain contributions from all spin S^^^. eigenstates. 
The hybridization TJhyb is treated as a perturbation, follow- 
ing Refs. l4lfl9ll24ll25l This leads to a master equation for the 
reduced density matrix pmoi in the Fock space of i/moi- 

The master equation still contains off-diagonal components 
of Pmoh corresponding to superpositions of molecular eigen- 
states. However, in the presence of non-commuting Zeeman 
and anisotropy terms in the Hamiltonian, any two states differ 
in the spin expectation value (Stot), which leads to different 
long-range magnetic fields. Thus the unavoidable interaction 
between the molecule and many degrees of freedom in the en- 
vironment (e.g., electron spins) should impart superselection 
rules^ ensuring that the dephasing of superpositions is rapid. 
For states with different charge, the Coulomb field also leads 
to strong superselection rules. ^1 

It is thus sufficient to consider the rate equations for the 
probabilities P,„ = (pmoi)mm of molecular states \m), 

Pni ^ ^ {,Rn — >mPn Rrn — >nPm\ (2) 

The stationary-state probabihties are determined by 

-Pm ' = 0. The transition rates Rm^n are written as a sum 
over leads and spin directions, Rm-,n — Ylaa ^niUw with 

To 

+ f{e„-e,r,-es^V/2)\C:J^], (3) 

where sl — 1, sr = — 1, I/tq = 27r |tp Dvuc/fi is the typical 
transition rate in terms of the density of states D (for one spin 
direction) of the leads and their unit-cell volume Vuc, V is 
the bias voltage, £,„ is the energy of state \m), f{x) is the 
Fermi function, and = (m|ccr|n) are transition matrix 
elements. The current through lead a can be expressed in 
terms of the rates and probabilities, 

^o. ^ -^a ^ ^ (j^m ^n) Rn — 'm (4) 

mncT 

where n,„ is the number of electrons for state |m). We are 
interested in the stationary state for which the current through 
both leads is equal. Therefore, we drop the subscript a and 
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insert the stationary-state probabilities P,!"^ . We will mainly 
analyze the differential conductance g = dl/dV. 

In principle, the solution of the equations Pm ^ = is sim- 
ple: The probabilities Pm form an eigenvector of a matrix 



A to zero eigenvalue, where A„ 



Rn 



for m n 



and A„ 



However, numerical diago- 
nalization often fails because the components of A, the rates, 
vary over many orders of magnitude due to the Fermi fac- 
tors in Eq. (|3]l. The ratio of very small rates can have a 
large effect on the probabilities, in particular for anisotropic 
magnetic molecules Ai2 Consequently, truncation errors in the 
very small rates can lead to large errors in the probabilities. 

An approach that avoids this problem would be highly wel- 
come. One such approach relies on the enumeration of all 
tree graphs on the network of molecular states connected by 
allowed transitions. =22i2iM However, this approach is not fea- 
sible for large molecular Fock spaces, since it requires sum- 
mation over of the order of N^~'^ tree graphs^ for each P,n ' , 
where N is the dimension of the space. 

The solution at T = is simpler and less susceptible to 
truncation errors because the Fermi factors are all either zero 
or unity. Thus there are no exponentially small rates, unless 
some matrix elements C,'^„ are exponentially small. This is 
not the case for generic angles 9 between magnetic field and 
easy axis. We here employ an approximation scheme to find 
the Pm ■* and the current /. The scheme relies on solving the 
rate equations and calculating the current exactly at T = and 
introducing broadening of the steps in P™ ^ and / to take finite 
temperatures into account. Details are presented in App. |A] 
The approximation is excellent at sufficiently low tempera- 
tures. It speeds up the calculation of g for 500 x 500 values 
of Vg and V and spin 5 = 2 by a factor of about 350. 



III. RESULTS AND DISCUSSION 

In the following, we discuss results for molecules with local 
spin 5 = 2 and 5=10. The smaller spin allows to exhibit the 
physics more clearly, since the number of relevant transitions 
is smaller As noted above, there are 2(25 + 1)^ transitions 
between states with zero and one electron, which gives 50 for 
5 = 2 and 882 for 5=10. These are the maximum possible 
numbers of differential-conductance peaks. For 5 = 2, we 
choose parameter values that allow to discuss the effects of 
interest but do not correspond to a specific molecule. 

For 5 = 10 we use realistic parameters for Mni2 calcu- 
lated by Park and Pederson^i employing density functional 
theory in the generalized-gradient approximation. We take 
A'2 = 0.0465 mcV, k = -0.00862 meV (taken from results 
for potassium doping), and J = 3.92 meV (from the energy 
difference of states with m = 21/2 and m = 19/2, respec- 
tively, where m is the quantum number of S^^^). K2 is close 
to the experimental value K2 = 0.056 mcV from Ref.lH 

The onside energy eq cannot be inferred from Ref. |27[ The 
LUMO-HOMO gap is Eg = 438 mcV.^^ Since Mnia is 
found to be more easily negatively doped, we can assume eo 
to lie closer to the LUMO. In experiments with a gate, the on- 




FIG. 1: (Color online) Differential conductance g — dl /dV as a 
function of gate voltage Vg and bias voltage V . Bright (dark) colors 
denote g > Q {g < Q). The parameters are 5 = 2, eo = 0, C/ = 10, 
H^^ H, = 0.05, K2 = 0.04, K = 0, J = 0.1, and T = 0.002 (all 
energies in meV). (a) Results obtained by exact solution of the rate 
equations, (b) Results of the approximation of App.lAl 



site energy is shifted to eq — eVg in any case. For single Mni2 
molecules, we will mainly consider the parameter region close 
to the crossing point where states with n = and n ~ 1 be- 
come degenerate. Finally, we choose a very large value for U 
so that double occupation of the LUMO is forbidden. 



A. Gate-voltage scans 

We start by presenting a typical differential-conductance 
plot as a function of gate voltage Vg and bias voltage V for 
spin 5 = 2 in Fig.[T] The magnetic field is applied at an an- 
gle = 45° relative to the easy axis. The comparison of the 
exact solution of the stationary-state rate equations with the 
approximation of App.lAl shows excellent agreement. 

The blue (medium gray) regions to the left and right of the 
crossing point in Fig.[T]correspond to Coulomb blockade (CB) 
with small current and electron numbers n = and n = 1, 
respectively. The plot would be mirror symmetric for V < 
0. The CB regions are delimited by strong peaks at the CB 
threshold, where the current increases rapidly to a large value 
of order c/tq. In the absence of internal degrees of freedom 
of the molecule this would be the only structure in the plots. 
However, the interaction of the tunneling electron with the 
local spin leads to inelastic tunneling processes, which cause 
the peaks at the CB threshold to split. 

The resulting fine structure is much more complex than for 
magnetic molecules without anisotropyi or in the absence of 
a strong magnetic fieldrii^ since there are many more allowed 
transitions in the present model due to non-commuting Zee- 
man and anisotropy terms. Each peak in g coiTesponds to one 
or more allowed transitions becoming energetically possible. 
The peak at the CB threshold is much stronger (the current 
step is much higher) than the peaks at higher bias voltage, 
since many additional transitions become available at the CB 
threshold, as illustrated by Fig.|2] 

Figure [T] also shows strong asymmetry between the fine 
structure on both sides of the crossing point. This is due to 
the multiplet of states with n = 1 being broader in energy 



4 




(a) n = n = l (b) n = n = I 



FIG. 2: Sketch of molecular energy multiplets (a) to the left of the 
crossing point in Fig.[T](ground state with n = electrons) and (b) to 
the right of the crossing point (ground state with n = 1). A few rep- 
resentative transitions are shown as double-headed arrows. Gener- 
ically, all transitions between states with n = and n — 1 have 
nonzero matrix elements. The bias voltage V is at the CB threshold 
in both cases. States that become populated at the CB threshold at 
T = are cross-hatched. In (a) some states with n = 1 do not 
become active at the CB threshold (white rectangle). 

than the one with n = 0. Figure |2] shows that for a ground 
state with n = 0, not all states with n = 1 become active 
(obtain nonzero probability at T = 0) at the CB threshold, 
whereas for a ground state with n ~ 1, all state with n = 
become active. Thus for the first case, to the left of the cross- 
ing point, there are more g peaks where additional states with 
n = 1 become active. 

Why then are there any additional lines to the right of the 
crossing point? Here, all n = states become active at the CB 
threshold. However, the probabilities and the current show 
steps when additional transitions become energetically pos- 
sible even if the final state of these transitions was already 
active. This mechanism leads to g peaks on both sides. 

Next, we consider parameters for Mni2. Figure [3] shows 
g{Vg, V) for a strong transverse field (Hz = 0) close to the 
crossing point between CB with n = and n = 1. Note that 
£0 has been set to zero, the zero of Vg is thus arbitrary. Figure 
|4{a) shows a sketch of special features in Fig. [3] 

The plots show two clear energy scales. For the strong mag- 
netic field (20 T) in Fig. [3] the Zeeman energy dominates over 
the anisotropy energy. The molecular states are nearly eigen- 
states of S'^Qj. In the following we use Fig. [3] as an example 
for the analysis of differential-conductance plots for complex 
molecules. We concentate on the case when the ground state 
has n = electrons, to the left of the crossing point. 

To facilitate the analysis. Fig. |4jb) shows the energy of 
low-lying molecular states vs. the expectation value (S'^^t), 
which is nearly a good quantum number Figure Hfc) shows, 
for each differential-conductance peak (current step) occur- 
ing at low bias voltage, the spin expectation value {S'^^^) for 
the initial (cross) and final (circle) state vs. the bias voltage. 
Only transitions starting from states with n = are included, 
these have negative slope in Fig. [3] Several transitions are 
marked with the same letters as in Fig. Ufa). These plots as- 
sume £0 ^ eVg = 25 mcV (the left edge of Fig. O, but the 
conclusions hold in general far left of the crossing point. 

The CB threshold corresponds to transition A in Figs.|4{a) 
and (c). Since all matrix elements between states with n = 
and n — 1 are nonzero, nearly all states belonging to the 
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FIG. 3: (Color online) Differential conductance g as a function of 
gate voltage Vg and bias voltage V for Mni2 in a magnetic field 
perpendicular to the easy axis. Bright (dark) colors denote g > 
(g < 0). The parameters are S = 10, e = 0, C/ = 10000, 
H:, = 2.32 (corresponding to 20 T), = 0, A"2 = 0.0465, 
K = -0.00862, J = 3.92, and T = 0.00863 (O.IK), where all 
energies are in meV. The approximation of App.lAlhas been used. 



low-energy ladder in Fig. Ub) assume nonzero probabilities 
even at T = as soon as transition A becomes energet- 
ically possiblei^ As noted above, this makes the threshold 
peak anomalously strong. 

The second strong line corresponds to transition B. With- 
out anisotropy, A and B would be the only visible transitions 
from n = to 71 = 1, since transition A (B) would coiTe- 
spond to a change in (S'fot) by +1/2 (—1/2) and these would 
be the only allowed changes. Also, the ladders of rt = 
and n = 1 levels would be exactly parallel so that all al- 
lowed transitions would have one of these two energies. In 
our case, the anisotropy leads to additional allowed transitions 
C, D,. . . etc. corresponding to changes of (S'toj) by approxi- 
mately —3/2,-5/2,..., respectively. 

The peaks show additional fine structure with a smaller en- 
ergy scale coming from the anisotropy. The CB-threshold 
peak is accompanied by additional peaks at slightly higher 
bias. The first visible one corresponds to transition Ai in 
Fig. Etc). Peaks B, C, D are accompanied by series of peaks 
at lower bias, reaching down to B', C, D'. These peaks coiTe- 
spond to arc-shaped series of transitions with approximately 
the same change in (S'fot). There are a few additional transi- 
tions with large changes in {Sf^^} in Fig.|4|c). These are not 
visible in Fig. [3] due to very small transition matrix elements. 

Several peaks, in particular B', show negative differential 
conductance (NDC). The origin is the following:- The current 
equals the charge e times a typical tunneling rate. This tunnel- 
ing rate is a weighted average over the rates of energetically 
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FIG. 4: (Color online) (a) g peaks corresponding to special transi- 
tions in Fig. [3] see text, (b) Molecular energy levels vs. spin expecta- 
tion value {Stot}, for — eVg = 25 mV. Open (solid) circles corre- 
spond to states with n — {n — 1) electrons in the LUMO. (c) Spin 
expectation values (Stot) of initial (crosses) and final (circles) states 
vs. bias voltage for the lowest-lying g peaks, for eo — eVg = 25 meV. 
Only transitions starting from states with ti = are shown. Special 
transitions are labeled as in (a). 



FIG. 5: (Color online) (a) Differential conductance g as in Fig. [3] 
except for a smaller magnetic field, Hx ~ 0.116 me V (IT). The 
approximation of App.|A]has been used, (b) Molecular energy levels 
vs. spin expectation value {Siot)^ for ^o — eVg = 19.5 mV. Open 
(solid) circles correspond to states with ?i = (n = 1). (c) Spin ex- 
pectation values (S'fot) of initial (crosses) and final (circles) states vs. 
bias voltage for the lowest-lying g peaks, for eo — eVg — 19.5 meV. 
Only transitions starting from states with n = are shown. 



possible transitions. But transition B' has a small rate due to 
small C,'^„, as numerical calculation shows. Therefore, the 
typical tunneling rate and the current decrease at B'. 

For comparison. Fig. |5ja) shows the differential conduc- 
tance for a smaller transverse field leading to comparable Zee- 
man and anisotropy energies. In this case, the molecular states 
are no longer approximate eigenstates of Sf^^, cf. Fig. HJb). 
This leads to a more complicated fine structure. It is still 
possible to attribute peaks to specific molecular transitions, as 
comparison with the map of transitions in Fig.|5jc) shows, but 
the spectrum and the transition energies look more random. 

Note that the fine structure splitting is now broader where 
the ground state has n ~ 1 electrons, opposite to the S = 2 
molecule, cf. Fig. [T] This stems from the smaller anisotropy 
A'2 + K < K2 in the n = 1 case, which makes the multiplet 
of accessible n = 1 states narrower in energy than the n = 
multiplet, as seen in Fig.|5jb)<^ 

B. Magnetic-field scans 

Our main topic is the interplay of magnetic field and 
anisotropy. Since the Zeeman and anisotropy terms in the 



Hamiltonian do not commute and are assumed to be compara- 
ble in magnitude, we expect the differential conductance g to 
depend significantly on the field. The natural way to study this 
is to plot g as a function of quantities characterizing the mag- 
netic field and possibly of bias voltage. Such magnetic-field 
scans could, in principle, also be done experimentally. 

Since a gate electrode is not required, magnetic-field scans 
could also be taken for monolayers of molecules sandwiched 
between conducting electrodes - or with an STM for single 
molecules or monolayers. For the theory to be applicable, the 
tunneling between molecule and electrodes must be made suf- 
ficiently weak to justify the sequential-tunneling approxima- 
tion, e.g., by molecular spacer groups or a thin oxide layeri^ 
The present subsection thus applies to single molecules and to 
monolayers of molecules with aligned easy axes. 

Figure|6ta) shows g as a function of magnetic-field compo- 
nents Hx and Hz for the case of spin S = 2 for fixed gate and 
bias voltages, eo — eVg = 0.2 meV and V = 0.8 mV, respec- 
tively. This plot shows g if one fixes Vg and V in Fig. [T] and 
varies the magnetic field. A peak in g appears whenever an al- 
lowed and energetically possible transition crosses the energy 
eV/2. Figure |6l a) can be extended to all values of H using 
the rotational symmetry of g around the z (easy) axis and the 
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FIG. 6: (Color online) (a) Differential conductance g as a function 
of the magnetic field components Hx (perpendicular to easy axis) 
and Hz (parallel to easy axis). Bright (dark) colors denote .g > 
{g < 0). The parameters are 5 = 2, eo - eVg = 0.2, U = 10, 
K2 = 0.04, K = 0, J = 0.1, T = 0.002 (all energies in meV), and 
V = 0.8 mV. The rate equations have been solved exactly, (b) <; as a 
function of the angle 6 between magnetic field and easy axis and bias 
voltage V for the S — 2 model in a magnetic field |H| = 0.1 meV. 
The other parameters are as in Fig.[T] The approximation of App.lAl 
has been used. 



reflection symmetry in the xy plane. 

Clearly, the symmetry of g(H) allows to determine the ori- 
entation of the molecule(s) from a transport measurement in 
a magnetic field. This is useful since the orientation is typi- 
cally poorly controlled and is not easy to determine in break- 
junction and electro migration experiments. 

The structure in Fig.|6ja) is richer than that in Fig.[T] Curves 
of large g are not straight and they vary in intensity (height of 
the current step) and can even vanish. They are not straight be- 
cause they correspond to differences of two eigenenergies of 
i/moi and the eigenenergies are complicated functions of H, 
since the Zeeman term does not commute with the anisotropy 
term. Nevertheless it is remarkable that the relatively simple 
model with spin S — 2 generates this complex structure. 

The g peaks in Fig.|6la) change in intensity with magnetic 
field since the transition matrix elements change. In the 
limiting case of = 0, i.e., field parallel to the easy axis, the 
Zeeman and anisotropy terms do commute. In this case S^^^ is 
conserved, and transitions changing S^^^ by values other than 
±1/2 are forbidden. Curves belonging to these transitions 
vanish for 0. 

Figure |6lb) shows g as a function of bias voltage V and the 
angle 6 between field and easy axis. Note again that several 
peaks vanish for 6^0, where the transitions become forbid- 
den. The plot shows that the strong peak at the CB thresh- 
old also depends on the magnetic-field direction. It is thus 
possible to switch the molecule between CB, with very small 
cuiTent due to cotunnelingr^ and a state with large current, 
leading to a large and anisotropic magnetoresistance at low 
temperatures. 

Figure|2ta) shows g{9, V) for Mni2. On first glance, sim- 
ilar features as in Fig. |6jb) are seen. However, an impor- 
tant new effect is at work here: Analysis of the 9 = case 
shows that the CB threshold for = corresponds to the 
transition from m = 10 to m = 21/2, where now rn is the 
quantum number of S^^^. The corresponding bias voltage is 
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FIG. 7: (Color online) (a) Differential conductance g as a function of 
the angle Q between magnetic field and easy axis and bias voltage V 
for Mni2 in a magnetic field |H| = 0.116 meV (1 T) at temperature 
T = 0.00863 meV (0.1 K). The other parameters are as in Figs. [3] 
and|5] The approximation of App.|A]has been used. The solid (green) 
curve denotes the CB threshold, which is invisible in t; for 6 > 
due to the presence of a blocking state, (b) Molecular energy levels 
vs. spin expectation value (Sfot) for Q = 0.3. Open (solid) circles 
correspond to states with n = Q (n = Y) electrons. The blocking 
state has n = and {Stot} ~ —10. Transitions discussed in the text 
are marked by arrows (most allowed transitions are not marked). 



VcB — 1-41 mV. However, Fig.lTJa) does not show a strong 
peak at this bias for small 6 > 0. Rather, the threshold ap- 
pears to be at about V = 1.64 mV. 

The molecular levels shown in Fig.|7Ib) help to understand 
this situation: the first possible transition at T = is the 
one from (S'tot) ~ 10 to {S^^^^) « 21/2, denoted by a solid 
curve in Fig. (Tja) and a solid aiTow in (b). For 9 = this 
leads to a large current associated with repeated transitions 
between these two states, which are the only allowed ones un- 
til the transition to (S^^^) « 19/2 becomes possible at higher 
bias. However, for 9 > there are nonzero matrix elements 
for transitions that decrease (S'tot) ^y about 3/2 (dashed ar- 
rows). Consequently, the molecule can reach the states with 
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(S'j^Qj) « 0. From there, it can relax to states with negative 
through transitions with large rates, in particular to the 
state with (S^^^) « —10. But this is essentially a blocking 
state, since the transitions to {S^^^) « —21/2 and —19/2 are 
energetically impossible and the only transitions that are ener- 
getically possible correspond to very large changes of (S'fot) 
(dotted arrows) and have extremely small matrix elements. 
This is thus an example of spin blockade. Note that cotim- 
neling transitions out of the blocking state to (S^^^) « —9 are 
possible, since the full potential difference eV is available for 
excitations. This leads to a very small current dominated by 
cotunneling, which would still be invisible in Fig.|7ja). 

To summarize, while transitions from positive to negative 
(S'^^Qj) happen only with a small rate, transitions in the op- 
posite direction occur with a still much smaller rate. Conse- 
quently, the molecule is nearly always in the (S^^^) w — 10 
state and the current is extremely small.— Thus there is no 
visible differential-conductance peak at the CB threshold. 

We conclude with two remarks: (i) The discontinuity of g 
at 6* = is irrelevant in practice, since it is impossible to per- 
fectly align the molecular easy axis with the field, (ii) Figure 
|7ja) shows that the CB threshold coincides with the strong g 
peak for 9 = it/2. In this case, Fig.|7jb) would be symmetric 
under spin inversion so that the states with (S^^^) « ±10 be- 
come degenerate ground states and there is no blocking state. 
This case was studied in Figs. |3] and |4] 
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FIG. 8: (Color online) Differential conductance per molecule av- 
eraged over all magnetic-field directions, 'g, as a function of gate 
voltage Vg and bias voltage V. Bright (dark) colors denote g > 
(g < 0). The parameters are identical to Fig. [T] except that the 
magnetic-field direction is not held fixed but is averaged over, and 
the same color scheme is used. 



C. Orientational disorder 

So far we have considered tunneling through single mole- 
cules or through monolayers of molecules with parallel easy 
axes.— Depending on the ligands, Mni2 can be nearly spheri- 
cal. In that case the distribution of orientations in a monolayer 
will be closer to random. 

To describe tunneling through a monolayer of randomly 
oriented magnetic molecules, we have to average the current 
and differential conductance over all orientations. Equiva- 
lently, we here average over all magnetic-field directions rel- 
ative to the molecular easy axis. One might expect that this 
smears out most of the fine structure. As we shall see, this 
is not the case. The averaged differential conductance g as a 
function of bias voltage V is rather complex due to van Hove 
singularities coming from extrema and crossings in the tran- 
sition energies as functions of magnetic field direction. 

The angle-averaged differential conductance per molecule 
is ^ = (l/47r) J dO dcf) (sin 6*) g, where 9, (p are the polar an- 
gles of the field, g shows peaks at bias voltages V correspond- 
ing to molecular-transition energies that depend on the angles. 
Compared to the conventional van Hove singularities of band 
theory, 5 corresponds to the density of states, V to the energy, 
and {9, 4>) to the wave vector. If the transition energies de- 
pended on both 9 and <f>, they would form a two-dimensional 
band structure in a "Brillouin zone" that is the surface of a 
sphere. Extrema in the transition energies would then lead 
to typical two-dimensional, step-like van Hove singularities. 
However, in our model the transition energies are indepen- 
dent of (j) so that 5 — (1/2) / dd (sin 6*) g. We obtain a one- 



dimensional band structure with an additional weight factor 
sin 6*. 

For an extremum at < 6* < tt, this leads to a one- 
sided singularity of the form 1 / y/\V — Vc\, typical for one- 
dimensional systems. If the extremum is at 6* = (and 9 ~ tt) 
and the transition is allowed there, the factor sin 9 reduces the 
singularity to a step. If the transition is forbidden for 9 = 0, 
the peak height in g is found to vanish as 9^, which leads to 
an even weaker singularity linear in |y — There are other 
cases of van Hove singularities, which are not present in nor- 
mal band structures and will be discussed below. 

It is the calculation of angle-averaged quantities for which 
the approximation scheme of App. lAl becomes crucial to re- 
duce the computational effort. All averaged quantities are cal- 
culated with this method. To approximate 

g = l f d9{sm9)g= f dug{u,V), (5) 
^ Jo Jo 

where u = cos 9, we first calculate g for fixed u at T = 0. We 
obtain a set of bias voltages Vi (u) and delta-function weights 
gi{u), which depend on u. We replace the integral by a sum 
over typically N = 5000 terms, 

1 ^ 

n— 1 i 

To obtain approximate results for T > we broaden the delta 
functions (or current steps) as described in App.lAl 

For orientation, we first show ^(Vg, V) for S* = 2 in Fig.|8] 
The plot should be compared to Fig. [T]for the same parameters 
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FIG. 9: (a) Differential conductance per molecule averaged over all 
magnetic-field directions, g, as a function of bias voltage V at three 
temperatures, 'g is plotted on a logarithmic scale. The parameters are 
as in Fig. [6{b). The inset shows g in the vicinity of the CB thresh- 
old on a linear scale. Several van Hove singularities are marked, (b) 
Upper part: bias voltage of g peaks as a function of magnetic-field 
angle 9, for the same parameters. Lower part: close-up of the vicin- 
ity of the CB threshold. The points corresponding to the van Hove 
singularities in (a) are marked by the same letters. 



except for fixed 9 = 45° in Fig.[T] Figure [8] shows that a lot 
of the fine structure survives the angular averaging. Plots of 
this type could not be obtained experimentally, due to the lack 
of a gate electrode for a monolayer. What can be measured 
is the differential conductance ^ as a function of bias voltage 
for given onsite potential eg. Large signals are expected if 
the particular molecule allows to reach the CB threshold with 
accessible bias voltages V. 

Figure ^a) shows 'g{V) for a monolayer of spin 5 = 2 
molecules in a magnetic field of magnitude |H| = 0.1 mcV. 
This plot corresponds to Fig. |6lb) averaged over 6 with the 
weight factor sin 9. We again see that the averaged differential 
conductance retains a lot of structure, in particular at very low 
temperatures. In the following, we analyze this structure in 
terms of van Hove singularities. 

To help with this analysis, Fig.|9jb) shows the bias voltages 
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FIG. 10: Differential conductance per molecule averaged over all 
magnetic-field directions, 'g, as a function of bias voltage V for a 
monolayer of Mni2 molecules with randomly oriented easy axes at 
two temperatures, g is plotted on a logarithmic scale. The inset 
shows ^ on a linear scale to exhibit the NDC region. Several van 
Hove singularities are marked, see text. 



of g peaks (current steps) as a function of magnetic-field angle 
9. These are the "bands" that lead to van Hove singularities. 
This is simply a map of g peaks occuring in Fig. |6jb). The 
first, strong peak in Fig.|9ja) corresponds to the CB threshold. 
Figures |6tb) and |9jb) show that the CB threshold is angle- 
dependent. Singularity A stems from the minimum of the 
transition energy and B from the maximum at 6* = 7r/2. These 
are both of the typical one-dimensional form 1 / y^\V — Vc\. 
The maximum at 6* = leads to singularity D, which is only a 
step, due to the weight factor sin 9. As noted above, the pres- 
ence of a step shows that the transition is allowed for 9 = 0. 
The 1 / ^y\V — Vc\ singularity E stems from the maximum of 
a transition not at the CB threshold. 

In addition, there is weaker step C, which cannot be at- 
tributed to an extremum at 6 = 0. Figure |9|b) shows that 
a transition can suddenly vanish when it intersects another 
one. This is a typical property of inelastic sequential tunneling 
through molecules, which is due to the initial state of one of 
the transitions being populated (at T = 0) only on one side of 
the crossing. One can show that the sum of weights (current 
step heights) of the transitions is analytic through the crossing, 
as a function of H^, H^, 9, or This means that a three- 
way crossing can be viewed as the superposition of a band 
for which the energy and the weight are analytic functions of 
9 and another band for which the weight is analytic but the 
energy shows a kink (change of slope). This kink leads to a 
type of van Hove singularity not present in band structures for 
weakly interacting electrons in solids. It shows up as a kink in 
the current and thus as a step in Singularity C in Fig.|9{a) 
is of this type, coming from the three-way crossing marked 
in Fig.|9tb). If both transitions have nonzero weight on both 
sides, part of the weight is also typically transferred from one 
band to the other, leading to the same type of singularity. 

Finally, we turn to a monolayer of Mni2 molecules {S = 
10) with random easy axes. We consider the case correspond- 
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ing to averaging g in Fig. |7ta) over angles. The result is the 
differential conductance plotted in Fig. [10] The most promi- 
nent features are the step-like singularities A and B. B stems 
from the maximum of the transition energy of the strong low- 
bias peak in Fig. |7ja) at about V ~ 1.64 mV. It is a step 
since the maximum occurs at 6* = 0. However, singularity 
A is also a step, although it clearly results from the mini- 
mum of that transition energy at 6* = 7r/2 so that we expect 
a pole, 1 / -^Z I y — T4 1 . This can be understood by extending 
Fig- 13 a) to angles 6* > 7r/2 using reflection symmetry. The 
minimum is in fact a crossing between the transitions from 
i^tot) = ±10 to ±21/2, i.e., the transitions out of the respec- 
tive ground state and blocking state. Therefore singularity A 
is not of the form of a band extremum but of a band crossing, 
i.e., a step. The anomalous exponent of singularity A is thus 
closely related to the spin blockade discussed earlier. 

The sharp peak X in Fig. [TO] is an artifact coming from 
molecules with 6* = 0, for which the CB threshold formally 
occurs at Vcb = 1.41 mV. Finally, the extended region of 
NDC in Fig. |7ja) leads to NDC even in the angle-averaged 
differential conductance, as seen in the inset of Fig.fTOl 



IV. SUMMARY AND CONCLUSIONS 

In summary, we have studied inelastic electron tunneling 
through molecules with a local magnetic moment and large 
uniaxial anisotropy in a strong magnetic field. Since the orien- 
tation of the molecules is often not well controlled in tunnel- 
ing experiments, we consider arbitrary angles between easy 
axis and field. Then the anisotropy and Zeeman terms in the 
Hamiltonian do not commute so that no component of the 
molecular spin is conserved. This lifts all spin selection rules 
for electron tunneling, leading to large numbers of allowed 
molecular transitions and consequently to many peaks in the 
differential conductance g ~ dl/dV. The resulting complex 
fine structure of Coulomb-diamond plots is already apparent 
for the relatively simple case of a local spin 5* = 2. 

As a concrete example, Mni2 molecules are studied. The 
large spin 5 = 10 leads to even more complex fine struc- 
ture. However, one can still attribute differential-conductance 
peaks to individual molecular transitions. It should be possi- 
ble to analyze experimental results in terms of specific molec- 
ular transitions ;/ one has a model for guidance. Even with- 
out detailed attribution of observed peaks, measurement of g 
in magnetic fields of various directions should allow to de- 
termine the orientation of the molecule relative to the leads, 
which is not directly accessable in break-junction or electro- 
migration experiments. 

We have considered three cases: Single molecules, molec- 
ular monolayers with aligned easy axes, and molecular mono- 
layers with random easy axes. For monolayers one does not 
have the advantage of a gate electrode. However, one can ex- 
tract similarly detailed information by varying the magnetic 
field. For randomly oriented molecules we have found that 
the averaging of transition energies over orientations leads to 
van Hove singularities in g. Besides the normal singularities 
from extrema of "bands," a novel type arises from crossings 



of transition energies. Analysis of the singularities can give 
rather detailed information on allowed vs. forbidden transi- 
tions in the limit of the easy axis aligned with the magnetic 
field and on the presence of blocking states (spin blockade). 

Detailed calculations of g for Mni2 with its many molec- 
ular states and in particular for monolayers with random ori- 
entation are made feasible by an approximation scheme for 
solving the rate equations. As further results, we predict large 
and highly anisotropic magnetoresistance at low temperatures, 
if the bias voltage is tuned close to the CB threshold, and neg- 
ative differential conductance, which survives even for ran- 
domly oriented monolayers of Mni2. 
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APPENDIX A: APPROXIMATE SOLUTION OF THE RATE 
EQUATIONS 

At T = we employ the following algorithm: For bias 
V ~ 0, the molecular ground states have probability 1/d, 
where d is the ground-state degeneracy, and the current van- 
ishes. For increasing V, the current remains zero in the 
sequential-tunneling approximation until eV/2 reaches the 
energy Ae of the first allowed transition starting from a ground 
state. Above this value, all states obtain nonzero probabilities 
that can be reached from a ground state directly or indirectly 
by allowed and energetically possible transitions. The net- 
work of these states is constructed by testing all possible tran- 
sitions. Then the matrix A above the threshold eV/2 = Ae is 
obtained from Eq. ([3]), taking into account that the Fermi fac- 
tors are / = or 1. A is then diagonalized. This is faster and 
more robust than in the general case because (a) the dimen- 
sion of A is smaller since it only contains active states and 
(b) A does not contain exponentially small components — all 
components are either exactly zero or given by matrix ele- 
ments. The resulting current is calculated from Eq. (|4|i. The 
value of Ae and the change (step height) in the current are 
recorded. Then we go over to the next bias for which eV/2 
equals a transition energy, extend the network of states by all 
states that can now be reached, obtain A, diagonalize it, and 
calculate the new current. These steps are repeated. 

Note that we only diagonalize ffmoi once to obtain the tran- 
sition energies and matrix elements. Furthermore, we only di- 
agonalize A once for each transition energy, since the current 
at r = is constant between steps. The result is a list of bias 
voltages Vi and associated heights gt of current steps. The gi 
are also the weights of delta-function peaks in g ~ dl/dV. 

To obtain approximate results at finite, but low, tempera- 
tures, we broaden the current steps so that their width is of the 
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order of fc^T. We write the current as 

^^^T.9^[f[-^^)-f[-^^)\^ (Al) 
where i enumerates the current steps. This particular broad- 



ening function is chosen since Eq. ( lAll ) is exact (assuming 
sequential tunneling) for the simplest possible model, i.e., a 
single orbital for a spin-less fermion. The form of the broad- 
ened delta functions in g follows trivially. The approximation 
is good if the separation in eV/2 between current steps is large 
compared to fc^T. 



* Electronic address: lctinmi@ku.edul 

' J. Park, A. N. Pasupathy, J. I. Goldsmith, C. Chang, Y. Yaish, J. 
R. Petta, M. Rinkoski, J. P. Sethna, H. D. Abrufia, P. L. McEuen, 
and D. C. Ralph, Nature (London) 417, 722 (2002). 

^ W. Liang, M. P. Shores, M. Bockrath, J. R. Long, and H. Park, 
Nature (London) 417, 725 (2002). 

^ J. Paaske and K. Flensberg, Phys. Rev. Lett. 94, 176801 (2005). 

* R Elste and C. Timm, Phys. Rev. B 71, 155403 (2005). 

^ H. B. Heersche, Z. de Groot, J. A. Folk, H. S. J. van der Zant, C. 
Romeike, M. R. Wegewijs, L. Zobbi, D. Bamca, E. Tondello, and 
A. Cornia, Phys. Rev. Lett. 96, 206801 (2006). 

* C. Romeike, M. R. Wegewijs, W. Hofstetter, and H. Schoeller, 
Phys. Rev. Lett. 96, 196601 (2006); Phys. Rev. Lett. 97, 206601 
(2006). 

^ C. Romeike, M. R. Wegewijs, and H. Schoeller, Phys. Rev. Lett. 

96, 196805 (2006). 
^ M.-H. Jo, J. E. Grose, K. Baheti, M. M. Deshmukh, J. J. Sokol, E. 

M. Rumberger, D. N. Hendrickson, J. R. Long, H. Park, and D. C. 

Ralph, Nano Lett. 6, 2014 (2006). 
' C. Timm and R Elste, Phys. Rev. B 73, 235304 (2006). 
'° R Elste and C. Timm, Phys. Rev. B 73, 235305 (2006). 
" M. N. Leuenberger and E. R. Mucciolo, Phys. Rev. Lett. 97, 

126601 (2006). 

A. Donarini, M. Grifoni, and K. Richter, Phys. Rev. Lett. 97, 
166801 (2006). 

M. Misiorny and J. Bamas, cond-mat/0610556 (unpublished). 

G. Gonzalez and M. N. Leuenberger, cond-mat/0610653 (unpub- 
lished). 

P. Elste and C. Timm, Phys. Rev. B (submitted), cond- 
mat/0611108. 

C. Joachim, J. K. Gimzewski, and A. Aviram, Nature (London) 

408, 541 (2000). 
" W. Hameit, Phys. Rev. A 65, 032322 (2002). 

Y. Xue and M. A. Ratner, in Nanotechnology: Science and 

Computation, edited by J. Chen, N. Jonoska, and G. Rozenberg 

(Springer- Verlag, Berlin, 2006). 
" K. Blum, Density Matrix Theory and Applications (Plenum, New 

York, 1981). 

H. Schoeller and G. Schon, Phys. Rev. B 50, 18436 (1994). 

^' J. Konig, H. Schoeller, and G. Schon, Europhys. Lett. 31, 31 
(1995). 

22 M. Turek and K. A. Matveev, Phys. Rev. B 65, 115332 (2002). 



H. Bruus and K. Plensberg, Many-body Quantum Theory in Con- 
densed Matter Physics (Oxford University Press, Oxford, 2004). 
A. Mitra, I. Aleiner, and A. J. Millis, Phys. Rev. B 69, 245302 
(2004). 

J. Koch, P. von Oppen, Y. Oreg, and E. Sela, Phys. Rev. B 70, 
195107 (2004). 

J. Koch and P. von Oppen, Phys. Rev. Lett. 94, 206804 (2005). 

" K. Park and M. R. Pederson, Phys. Rev. B 70, 054414 (2004). 

~^ J. R. Friedman, M. P. Sarachik, J. Tejada, and R. Ziolo, Phys. Rev. 
Lett. 76, 3830 (1996). 

2' K. M. Mertes, Y. Suzuki, M. P. Sarachik, Y. Myasoedov, H. Shtrik- 
man, E. Zeldov, E. M. Rumberger, D. N. Hendrickson, and G. 
Christou, Solid State Commun. 127, 131 (2003). 

^" B. Fleury L. Catala, V. Hue, C. David, W. Z. Zhong, P Jegou, 
L. Baraton, S. Palacin, P. -A. Albouy, and T. Mallah, Chem. Com- 
mun. (Cambridge) 2005, 2020. 

^' W. H. Zurek, Phys. Rev. D 26, 1862 (1982). 

32 T. L. Hill, J. Theor. Biol. 10, 399 (1966). 

33 J. Schnakenberg, Rev. Mod. Phys. 48, 571 (1976). 

3" R. K. P Zia and B. Schmittmann, J. Phys. A: Math. Gen. 39, L407 
(2006). 

3^ The higher-energy multiplet of states with n — 1, which is split 
from the lower multiplet by an energy of order JS, also becomes 
occupied due to tiny matrix elements between n = states with 
large negative langleStot) and n = 1 states with large posi- 
tive langleStot) ■ However, the probabilities are negligible for the 
present discussion. 

3'' There is another multiplet of n = 1 states at higher energies, split 
from the low-energy multiplet by a large energy of the order of 
JS. These states are not populated for the bias voltages consid- 
ered here. 

3^ X. H. Qiu, G. V. Nazin, and W. Ho, Phys. Rev. Lett. 92, 206102 
(2004). 

3^ A significant number of electrons have to tunnel through the 
molecule for it to get from the ground state to the state with 
(S'tot) ~ ^10, but after that it is essentially trapped and the sta- 
tionary current is very small. 

3' This follows because the current at T = is analytic, as long as no 
transition is crossed, since the Fermi functions are then constant 
and the matrix elements are analytic functions of magnetic field 
and gate voltage. 



